


#### APPENDIX C5: interactions

# load, combine data
library(rio)
library(tidyverse)

# estimate models
library(brglm2)

# cluster standard errors
library(lmtest)
library(sandwich)

# generate tables
library(stargazer)

# generate plots
library(sjPlot)
library(scales)
library(ggpubr)




# LOAD DATA
glm_df <- import("glm_data.csv") %>%
  filter(iso3c!="TTO" & iso3c!="GUY" & iso3c!="SUR") %>%
  filter(!is.na(approval) & !is.na(opp1vote) & !is.na(number_protests))



## TABLE C8: INTERACTIONS
brglm_interact1 <- glm(doc ~ approval * opp1vote + number_protests + previous_doc + duck + partyage + polity2 + left_exec + any_election +
                      minister_is_technocrat_mainstream + resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + 
                      gdp_growth_lag + imf_program + crisis + as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
              family = binomial(link="logit"), 
              method = "brglmFit",
              data = glm_df)

brglm_interact1_r <- brglm_interact1 %>%
  coeftest(vcovHC(brglm_interact1, type = 'HC0', cluster =  ~ iso3c))


brglm_interact2 <- glm(doc ~ approval * number_protests + opp1vote + previous_doc + duck + partyage + polity2 + left_exec + any_election +
                       minister_is_technocrat_mainstream + resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                       as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
                       family = binomial(link="logit"), 
                       method = "brglmFit",
                       data = glm_df)

brglm_interact2_r <- brglm_interact2 %>%
  coeftest(vcovHC(brglm_interact2, type = 'HC0', cluster =  ~ iso3c))


brglm_interact3 <- glm(fund_doc ~ approval * opp1vote + number_protests + previous_fund_doc + duck + partyage + polity2 + left_exec + any_election +
                       minister_is_technocrat_mainstream + resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                       as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
                       family = binomial(link="logit"), 
                       method = "brglmFit",
                       data = glm_df)

brglm_interact3_r <- brglm_interact3 %>%
  coeftest(vcovHC(brglm_interact3, type = 'HC0', cluster =  ~ iso3c))


brglm_interact4 <- glm(fund_doc ~ approval * number_protests + opp1vote + previous_fund_doc + duck + partyage + polity2 + left_exec + any_election +
                       minister_is_technocrat_mainstream + resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                       as.factor(iso3c) + as.factor(quarter) + as.factor(year), 
                       family = binomial(link="logit"), 
                       method = "brglmFit",
                       data = glm_df)

brglm_interact4_r <- brglm_interact4 %>%
  coeftest(vcovHC(brglm_interact4, type = 'HC0', cluster =  ~ iso3c))


# generate tables
stargazer(brglm_interact1_r,brglm_interact2_r,brglm_interact3_r,brglm_interact4_r, #type="text",
          omit=c("iso3c","year","quarter"), no.space = T,
          covariate.labels = c("Executive approval (perc)","Opposition vote share (perc)","Number of protests",
                               "Previous document = 1","Previous fund document = 1","Lame duck = 1", "Party age",
                               "Democracy","Left executive = 1",
                               "Election quarter = 1","Mainstream minister = 1",
                               "Resource rents (perc)","Field discovery = 1","Oil price (USD, log)",
                               "GDP per capita (1,000 USD)","GDP growth (perc)","IMF program = 1", "Crisis = 1",
                               "Executive approval X Opposition vote share", "Executive approval X Number of protests"), 
          dep.var.labels = c("Any document","Fund document"))

# retrieve n, aic, ll, etc
stargazer(brglm_interact1,brglm_interact2,brglm_interact3,brglm_interact4, #type="text",
          omit=c("iso3c","year","quarter"), no.space = T,
          dep.var.labels = c("Any document","Fund document"))





## PLOT 1: EXECUTIVE APPROVAL X OPPOSITION VOTE SHARE - ALL DOCS
brglm_plot1 <- glm(doc ~ approval * opp1vote + number_protests + previous_doc + duck + partyage + polity2 + left_exec + any_election +
                         resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                     #as.factor(iso3c) + 
                     quarter + year, 
                       family = binomial(link="logit"), 
                       method = "brglmFit",
                       data = glm_df)

plot_appendix1 <- plot_model(brglm_plot1, type = "pred", terms = c("approval[all]", "opp1vote[0]")) + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = expression(paste("(a) ", italic("Opposition Vote Share"), " at Its Minimum (0%)"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,1), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix2 <- plot_model(brglm_plot1, type = "pred", terms = c("approval[all]", "opp1vote[21.02]")) + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = expression(paste("(b) ", italic("Opposition Vote Share"), " at Its Median (21.02%)"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,1), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix3 <- plot_model(brglm_plot1, type = "pred", terms = c("approval[all]", "opp1vote[22.85]")) + # mean value of opp1vote
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = expression(paste("(c) ", italic("Opposition Vote Share"), " at Its Mean (22.85%)"))) + scale_x_continuous(labels = label_percent(scale = 1))  + 
  scale_y_continuous(lim = c(0,1), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix4 <- plot_model(brglm_plot1, type = "pred", terms = c("approval[all]", "opp1vote[52.20]")) + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = expression(paste("(d) ", italic("Opposition Vote Share"), " at Its Maximum (52.20%)"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,1), labels = label_percent(scale = 100, accuracy = 1))

ggsave("interaction-plot1.pdf", ggarrange(plot_appendix1,plot_appendix2,plot_appendix3,plot_appendix4), width = 9.75, height = 8)


## PLOT 2: EXECUTIVE APPROVAL X NUMBER OF PROTESTS - ALL DOCS
brglm_plot2 <- glm(doc ~ approval * number_protests + opp1vote + previous_doc + duck + partyage + polity2 + left_exec + any_election +
                         resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                     #as.factor(iso3c) + 
                     quarter + year, 
                       family = binomial(link="logit"), 
                       method = "brglmFit",
                       data = glm_df)

plot_appendix5 <- plot_model(brglm_plot2, type = "pred", terms = c("approval[all]", "number_protests[0]")) + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = expression(paste("(a) ", italic("Number of Protests"), " = 0"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.4), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix6 <- plot_model(brglm_plot2, type = "pred", terms = c("approval[all]", "number_protests[1]"))  + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = expression(paste("(b) ", italic("Number of Protests"), " = 1"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.4), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix7 <- plot_model(brglm_plot2, type = "pred", terms = c("approval[all]", "number_protests[2]"))  + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = expression(paste("(c) ", italic("Number of Protests"), " = 2"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.4), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix8 <- plot_model(brglm_plot2, type = "pred", terms = c("approval[all]", "number_protests[3]"))  + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Any Document"))),
                         title = expression(paste("(d) ", italic("Number of Protests"), " = 3"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.4), labels = label_percent(scale = 100, accuracy = 1))


# ggsave("interaction-plot2.pdf", ggarrange(plot_appendix5,plot_appendix6,plot_appendix7,plot_appendix8), width = 9.75, height = 8)



## PLOT 3: EXECUTIVE APPROVAL X OPPOSITION VOTE SHARE --- FUNDS
brglm_plot3 <- glm(fund_doc ~ approval * opp1vote + number_protests + previous_fund_doc + duck + partyage + polity2 + left_exec + any_election +
                     resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                     #as.factor(iso3c) + 
                     quarter + year, 
                   family = binomial(link="logit"), 
                   method = "brglmFit",
                   data = glm_df)

plot_appendix9 <- plot_model(brglm_plot3, type = "pred", terms = c("approval[all]", "opp1vote[0]")) + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = expression(paste("(a) ", italic("Opposition Vote Share"), " at Its Minimum (0%)"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,1), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix10 <- plot_model(brglm_plot3, type = "pred", terms = c("approval[all]", "opp1vote[21.02]")) + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = expression(paste("(b) ", italic("Opposition Vote Share"), " at Its Median (21.02%)"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,1), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix11 <- plot_model(brglm_plot3, type = "pred", terms = c("approval[all]", "opp1vote[22.85]")) + # mean value of opp1vote
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = expression(paste("(c) ", italic("Opposition Vote Share"), " at Its Mean (22.85%)"))) + scale_x_continuous(labels = label_percent(scale = 1))  + 
  scale_y_continuous(lim = c(0,1), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix12 <- plot_model(brglm_plot3, type = "pred", terms = c("approval[all]", "opp1vote[52.20]")) + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = expression(paste("(d) ", italic("Opposition Vote Share"), " at Its Maximum (52.20%)"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,1), labels = label_percent(scale = 100, accuracy = 1))

ggsave("interaction-plot3.pdf", ggarrange(plot_appendix9,plot_appendix10,plot_appendix11,plot_appendix12), width = 9.75, height = 8)


## PLOT 2: EXECUTIVE APPROVAL X NUMBER OF PROTESTS --- RULES
brglm_plot4 <- glm(fund_doc ~ approval * number_protests + opp1vote + previous_fund_doc + duck + partyage + polity2 + left_exec + any_election +
                     resource_rents_lag + discovery + real_oil_price_log + gdp_pc_constant_lag + gdp_growth_lag + imf_program + crisis +
                     #as.factor(iso3c) + 
                     quarter + year, 
                   family = binomial(link="logit"), 
                   method = "brglmFit",
                   data = glm_df)

plot_appendix13 <- plot_model(brglm_plot4, type = "pred", terms = c("approval[all]", "number_protests[0]")) + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = expression(paste("(a) ", italic("Number of Protests"), " = 0"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.4), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix14 <- plot_model(brglm_plot4, type = "pred", terms = c("approval[all]", "number_protests[1]"))  + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = expression(paste("(b) ", italic("Number of Protests"), " = 1"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.4), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix15 <- plot_model(brglm_plot4, type = "pred", terms = c("approval[all]", "number_protests[2]"))  + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = expression(paste("(c) ", italic("Number of Protests"), " = 2"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.4), labels = label_percent(scale = 100, accuracy = 1))
plot_appendix16 <- plot_model(brglm_plot4, type = "pred", terms = c("approval[all]", "number_protests[3]"))  + 
  theme_minimal() + labs(x = "Executive Approval", y = expression(paste("Predicted Probability of ", italic("Fund Document"))),
                         title = expression(paste("(d) ", italic("Number of Protests"), " = 3"))) + scale_x_continuous(labels = label_percent(scale = 1)) + 
  scale_y_continuous(lim = c(0,0.4), labels = label_percent(scale = 100, accuracy = 1))


ggsave("interaction-plot4.pdf", ggarrange(plot_appendix13,plot_appendix14,plot_appendix15,plot_appendix16), width = 9.75, height = 8)
